The liquid-glass transition of silica 
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We studied the liquid-glass transition of Si02 by means of replica theory, utilizing 



' an effective pair potential which was proved to reproduce a few experimental features 
> ■ 

■ of silica. We found a finite critical temperature Tq, where the system undergoes a 

^ : 

, phase transition related to replica symmetry breaking, in a region where experiments 

00 

O I do not show any transition. The possible sources of this discrepancy are discussed. 



I. INTRODUCTION 



B 

■ According to the behavior of their transport coefficients, supercooled liquids were classi- 

O . fied in fragile and strong liquidsB. Close to the experimental glass temperature T^l, where the 
system is no longer able to equilibrate on experimental time scales, the viscosity n in strong 

X ■ 

;_j ■ liquids seems to be in a good agreement with the Arrhenius law exp(yl/T), while in the frag- 

ile ones it is more correctly described by the Vogel-Fulcher-Tamman law exp(y4/(T — Tq)). 
As the dramatic slowing down of dynamics prevents to achieve equilibrium below Tg, the 
critical temperature Tq can only be inferred by extrapolating at lower temperatures the avail- 
able equilibrium data. However, this fitting procedure is hardly conclusive, as the viscosity 
changes of 14 order of magnitude in a small window of temperature^. 

In such a classification silica {Si02) should belong to the class of strong liquid, though 
deviations from the pure Arrhenius plot above Tg have been recently observed!, supporting 
speculations about a fragile-strong transition at temperatures higher than Tg, that for silica 
is known to b J ~ 1U6K. 

The amorphous network observed in supercooled silica has been recently shown! to arise 
in molecular dynamics simulation when using the effective two body potential recently in- 



2 



troduced by van Best, Kramer and Van SantenQ (BKS). However, due to the hmitations of 
numerical simulations, the transport coefficients of BKS silica were studied only in a regime 
of temperature well above Tg, hence the equilibrium behavior at low temperature is still an 
open issue. Controversial results indeed arose when the available data were extrapolated at 
lower temperatures and the possibility that the BKS potential rather implies a fragile behav- 
ior cannot ruled out. As a matter of fact, fitting data with a VFT law yields Tq = 2525idi, 
while the power law fit 77 oc (T — TcY, motivated by the predictions of mode-coupling the- 
ory!, yields the dynamical temperature Tc ~ 33OOi^i'00. On the other hand, the Arrhenius 
behavior seems to be recovered taking into account only the few data available at the lowest 
temperatures, suggesting the possibility of a fragile-strong transition at T^. Furthermore, the 
study on the configurational entropy S^, shows that the Adam-Gibbs relation rj oc 1/TSc 
is verified over 4 decades and suggests that the fragile-strong transition might be related to 
a change of the energy landscape properties at low temperatures. 

Since the reliability of the BKS potential in describing the low temperature behavior 
of the viscosity of silica is still far from being established, we believe it is instructive to 
apply the replica approachSElllllli to this case, looking for a possible critical temperature 
To from a different route, namely as the point where the system undergoes a thermodynamic 
transitionlll related to a replica symmetry breaking^. 



II. CORRELATED LIQUID PHASE 

Let us assume that the glassy dynamics is the signature of the partitioning of the whole 
phase space in an exponentially large number A/'(/) of free energy minima {valleys in the 
following), whose exploration in the supercooled phase is still allowed due to hopping pro- 
cesses between them. The entropic contribution coming from the huge number of different 
valleys is called complexity: 

S = \ogX{f)/N (1) 

The Valleys' must be thought as somewhat more complicated objects^ than the basins of 
the 'inherent structures'liiS, as their number can be badly different^. 

Generally speaking, a valley is stable against both an infinitesimal displacement of an 
extensive number of particles and a finite rearrangement of a non-extensive number of them. 
It has been arguedil that the complexity in a finite dimensional system should be thought 
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as a time-dependent quantity. The thermodynamic approach0'Ei'Ei rather assumes that the 
valleys are metastable states with a macroscopic time-life, whose existence is limited to the 
region below T^, as suggested by the analogy with generalized (mean field) spin glasses!!. 
We shall exploit this assumption in order to make predictions about the supercooled phase, 
the glass temperature and the glass phase. 

The computation of the complexity T) in a given thermodynamic state of the system 
(hereafter represented by its temperature T) can be settled by using the method proposed 
by MonassonEl, which amounts to compute the free energy of a system composed by m 
identical clones (copies) of the original system. 

In the liquid phase, upon introducing a suitable coupling, two arbitrary clones can be 
forced to stay in the same valley [correlated liquid phase), but when the coupling vanishes 
they are always found to be in different valleys [uncorrelated liquid phase). On the other 
hand, in the glass phase, the clones remain in the same valley even after sending the coupling 
term to zer 00. Generally speaking, one needs to compute the free energy density of the 
cloned system in its correlated liquid phase. The information on the complexity of the 
cloned system is obtained when one computes the analytic continuation of the correlated 
liquid free-energy density 0(m) to real values of m, by means o 

T dm ^ ' ^ ^ ^ 

Given this very general framework, the aim is to compute the potential 0, depending both 
on the standard thermodynamic parameters (temperature, density, chemical potential, etc..) 
and on the parameter m, once the pair potential V{r) has been given. Three different ways 
to compute the free energy of a replicated system were introduced0'i'0, each relying on 
a different set of approximations. The results shown in the following were obtained within 
a harmonic approach, which only requires the pair correlation functions g{r) computed 
in the liquid phase. The position x^""^ of the a-th clone of the i-th particle is written as 
x^""^ = Zi + u\"'\ isolating the centers of mass Zi. By assuming that the displacements u^""^ 
are 'small' in the correlated liquid phase, the partition function (in d dimensions) is: 

= m^'^/^v^^'^^™-^^ Zu,{T/m) (exp (-^^ Tr logT^/r) )^^^ , (3) 

where Ziig{T/m) is the partition function of the liquid phase at the temperature T/m and 
7i is the Hessian matrix of the system. The average < . . . > is taken over the equilibrium 
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configurations of tlie variables {zi}. Introducing tlie quenclied approximation^: 



^exp (-^ Trlog7^/T))^^^ ~ exp (-^ (Tr log7^/^)^/„) (4) 
tlie free energy 0(m) is: 



(m) = fHq{T/m) + ^\{m- 1) {TilogH/T)^, - d (m - 1) In 27r - d log ml (5) 



2m 



The computation of (Tr logTi/T)^,^,^ is simplified by treating the fluctuations of the diagonal 
term in the Hessian matrix in a perturbative way0 which is correct in the high density limit, 
while the off-diagonal term can be taken into account non-perturbatively00'0. 
Let us briefly recall the scenario and the predictions of the theory: 

• In the liquid phase the physical quantities are obtained in the limit m ^ 1. In that 
case formula (^) simply reads: 

T.{T) = Su,{T)-Ssoi{T), (6) 

Siiq being the entropy of the liquid and Ssoi the entropy of a disordered harmonic solid. 

• If a temperature To where complexity vanishes actually exists, there the ergodicity is 
broken and the system undergoes the thermodynamic glass transition. 

• Choosing as a suitable order parameter of the phase transition the correlation function 
between the different clones, the transition turns out to be discontinuous. Surprisingly 
enough from the thermodynamic point of view it is a second order phase transition, 
the heat capacity jumping downwards from liquid to crystal-like values. 

• In the glass phase the eq. at m = 1 yields a negative complexity, which is mean- 
ingless. Since in the m — T plane one still finds out a critical line m*(T) (lesser than 
1) where S = 0, the replica approach takes 0(m*,T) as the physical free energy. 

Finally, it is interesting to note that the replica scenario shares many features with the old 
Adam-Gibbs-Di Marzio one0. 



III. RESULTS 



We analytically studied the glass transition of Si02, when the temperature is lowered 
and the density is kept fixed at the value p = 2.36(?/cm^, which is close to the experimental 
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density at ambient pressure0'0. We studied the correlated hquid phase within the harmonic 
approach sketched in the previous section, describing the interactions by means of the BKS 
potential: 

V-\r) = ^ + A,,e-^«' ^ - ^ + (7) 

where Q,,> = q,q,>e^ with qsi = 2.4, go = -1.2, = 1602.19/(47r8.8542) eV ■ A, the other 
parameter values being reported in Tab. (|I[). Since the potential has no lower bounds when 
r ^ 0, we followed the procedure of numerical simulations adding a short-ranged 24 — 6 
Lennard- Jones term, which does not affect the thermodynamic quantitiesii. 

In order to obtain the g{r) we resorted to standard theory of liquidsii. Among the 
others, the so-called Hypernetted Chain Approximation (HNC) turns out to be a rather 
reliable approach in describing fluids with electrostatic long range interactions such as silica. 
In the HNC free energy the infrared divergence due to Coulomb term is exactly canceled 
out by the resummation of the divergences of the virial expansion^. Unfortunately, the 
HNC closure is not powerful enough to study the correlated liquid phase. As a matter of 
fact, a thermodynamic inconsistency between the two different routes for computing the 
compressibility makes the integral 

Xi = - + P [d'rigir)-!) (8) 

to increase far more rapidly than the quantity 

when going to low temperatures. This prevented to study the whole correlated liquid phase, 
since it is not possible to obtain the g{r) in HNC approximation for temperatures below 
~ 3200ir. 

A very simple procedure to deal with this difficulty is to get the results at very low 
temperatures by extrapolating the ones obtained at the temperatures allowed within the 
HNC approach. Since such a procedure might appear not satisfactory enough, we have also 
chosen to utilize a different approximation interpolating between HNC and Mean Spherical 
Approximation (MSA) introduced by Zerah & Hansenll (ZH), which retains the advantages 
of the HNC approximation extending its reliability at much lower temperatures, in such a 
way that no extrapolation is required. It is interesting to note that the extrapolated HNC 
results turn out to be very similar to the ones obtained by means of the ZH approach, as 
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we shall show. Finally let us point out that, since the short range behavior is ruled by the 
MSA closure while the long range one is still described by the HNC, the cancellation of the 
infrared divergence described above is unaffected by the interpolation. 

More in detail, we considered the following integral equations for the g'^''^'{r): 

ec'p'gAr) = exp (-/?<)(r)) 1 + ^ ^'^ ^[J^^^' ^^^^ (10) 



where V(^)(r) = V^" (r) — V^^^ for r < rmin and zero otherwise, while V^^^{r) = V" (r) — 
V('^j(r). Moreover w^^'{r) = J dr'h^^'(r — r')c"'(r'), where c^'^'(r) is the direct correlation 
function and h'^^' [r] = g^^' {r) — 1. The interpolation is achieved by means of the function 
/(r) (here the same for all the components). In the regions where /(r) ~ 1 the eq. (p!oD 
gives the well-known HNC equations, while the MSA ones are recovered for /(r) ~ 0. As 
we made the choice /(r) = 1 — exp(— r/a) the interpolation depends only on the parameter 
a. Its value can be suitably chosen in order to strongly reduce the difference between the 
two compressibilities. 

We numerically solved the equations obtaining the pair correlation functions in two dif- 
ferent cases. In the former we solved the HNC equations, corresponding to the a — > limit 



in eq. (10), while in the latter a has been chosen in such a way to reduce the difference 
Xi — X2- More precisely, we let a to take the constant value 0.7, which has been checked to 
be the average value minimizing the difference between the two compressibilities down to 
~ 4000-/^. This solution of the integral equations in the a = 0.7 case gave our 5'ee'(r) in the 
ZH approximation. 

In order to solve the integral equations we discretized them, introducing a cut-off L in 
the real space and a mesh size a. More in detail, we are showing results obtained when 
the cut-off in the real space has the value L = 32 A, considering respectively 4096 points 
(HNC) and 2048 points (ZH), checking that they do not change significantly below Tc when 
L is doubled, though we still observe small scale oscillations on the correlation functions 
depending on the mesh size. 

In figs.(|l],0J^) we show the pair correlation obtained in the two different approximations, 
comparing them at two different temperatures {AOOOK, 27507^") with numerical data pre- 
viously published by Horbach and Let us note that the positions of maxima and 
minima are quite correctly reproduced but there is an error even larger than 10% which 
increases when lowering the temperature on the values of the first peak and of the first 
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FIG. 1: The analytical Si — Si pair correlation functions compared with numerical data by Horbach 
and Kob. 
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FIG. 2: The analytical O — O pair correlation functions compared with numerical data by Horbach 
and Kob. 

minimum, though they roughly cancel out when considering the integral. Moreover the ob- 
tained g{r) show the already mentioned small scale oscillations which are related to finite 
size effects in the momenta domain. Therefore the agreement with simulations seems to be 
not as good as in the previously considered cases 

As stated before, the HNC results at very low temperatures can only be inferred by 
extrapolating the high temperature ones. That procedure is simplified by the theoretical 
argument given by Rosenfeld & Tarazonaii yielding a T~^/^ behavior for the entropy of 
the liquid (apart from the thermal contribution). Interestingly, also in the ZH approach we 
found a Sug depending linearly on T~^/^ in the whole low temperatures region. On the other 
hand, the temperature dependence of the disordered solid entropy is well approximated by a 
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FIG. 3: The analytical Si — O pair correlation functions compared with numerical data by Horbach 
and Kob. 
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FIG. 4: The liquid and solid entropies as obtained by extrapolating HNC results with the functions 
(iHq + hiiqT~'^/^ and agoi + &so« log T respectively. Here the thermal contribution and the correct 
normalization constant (giving S" = when the phase space volume is equal to h^) are added to 
both of the entropies. 

logarithmic behavior, as expected when only vibrations are taken into account. In fig.(^) we 
show both the entropies obtained by the HNC solutions and the results of the extrapolation 
procedure. 

The main result of our work is presented in fig. (1), showing the complexity of the system 
in the liquid phase, computed both in the HNC and ZH approaches utilizing the eq.@ 
when 0(m) is given by the eq. @. As anticipated, the two curves gives roughly the same 
To ~ 21007^ yielding a rather sound prediction, which should not strongly depend on the 
particular way chosen to describe the liquid phase. 
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FIG. 5: The complexity as obtained within the HNC and the ZH approximation respectively. 
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FIG. 6: The behavior of m as function of the temperature in the glassy phase. 



In the glass phase we want to focus on the predicted behavior of the thermodynamics 
parameter m* as a function of the temperature, shown in fig.(^. As explained before, m* 
has been found by maximizing the free energy in the correlated liquid phase with respect 
to m. This thermodynamic parameter has a linear behavior, shared by all the glassy systems 
to our knowledge, ranging from fragile glasseJiiB to generalized spin glassesii, implying that 
the quantity T/m* remains nearly constant in the whole glassy phase. The analogy with 
generalized spin-glasses suggests that m* might be related to the relative height of the two- 
peaks equilibrium probability distribution of distances between replicas. Some numerical 
evidence for this quantity being non trivial was obtained!^. 
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IV. DISCUSSION 

The study of supercooled BKS silica, by means of the replica approach dealing with 
.™c..al g,... ..oduced ,y M..d and wMc. u.„.e. ..p,e ..eo. ,n 

order to obtain the pair correlation functions, yielded a liquid-glass phase transition at a 
temperature Tq ~ 2100/^. This result is quite odd, since it contradicts the experimental 
result implying no thermodynamical transition down to the experimental glass temperature 
Tg = 14467^', hence a few comments are strongly required. 

As it should be clear from the previous discussion on the behavior of viscosity in the 
BKS silica, our theoretical findings are not ruled out on numerical groundsES, as MD 
simulations have not explored the region where temperatures are much lower than the Mode- 
Coupling one. Because of that a nearly-Arrenhius behavior (covering no more than two 
decades) has been observed only taking into account a few low temperature data, whereas 
a previous extrapolation of datai was compatible with a divergence at a finite temperature 
To larger than Tg. Therefore one could argue that the too large critical temperature found 
in this thermodynamic approach is nothing but an artifact of the BKS potential which is 
not completely faithful to the actual low temperature behavior of real silica. 

On the other hand, as we have already stressed, the g{r) obtained from approximate 
theory of liquids show some difference from the ones obtained by MD simulations utilizing 
the BKS potential. We want to recall here that both the approximations involved (HNC and 
ZH) satisfy the T~^/^ behavior for the entropy of the liquid in the whole low temperature 
phase, while the claimed fragile-strong transitionEl at points to the failing of this law. 
This discrepancy could be the source of a severe overestimatation of the value of the actual 
critical temperature. 

Hence we are led to consider 21007^ rather as a rough estimate of the actual glass temper- 
ature. As a matter of fact, experiments cannot rule out the possibility of a thermodynamic 
transition at a temperature lower than Tg. 

Furthermore, it is worthwhile to note that the approximations involved in the computa- 
tion of correlated liquid free energy are quite crude. This could lead to an overestimation of 
the critical temperature as well. 

An intriguing issue is the possibility to gain an insight into the existence of a finite Tq 
by means of measurements of the violation of the fluctuation-dissipation relation in off- 
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equilibrium experiments. It has been suggested the existence in structural glasses of an off- 
equilibrium regime, called 'aging', reached by a system driven out-of-equilibrium by letting 
it to suddenly cross the transition point. 

It has been also argued that in the aging regime the fluctuation-dissipation relation could 
still be utilized to define an effective temperature T^ff which does not depend on the age of 
the system0'ii. On very general grounds, it has been moreover conjectured the equivalence 
for disordered systems between the ratio T/Tgff and the parameter describing the replica 
symmetry breaking, that at one step level is nothing but the function m*(T) computed in the 
thermodynamic approach^. Of course, the existence of Tg makes the experimental situation 
much more complicated, as for practical purposes the system falls out-of-equilibrium when 
crossing Tg. At the moment, it is unclear how to link an effective temperature measured 
between Tg and Tq to some purely thermodynamic computation. 

Nevertheless, a few attempts on different systems have been already performedBS or are 
in progressii, and we believe that even on silica those kind of experiments could be useful 
to clarify many open points. 
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TABLE I: Parameters of the BKS potential 





A (<A0 


B (A-i) 


c (cv ■ A-«) 


SiSi 


0.0 


0.0 


0.0 


00 


4.87318 


2.76 


175.0 


SiO 


18003.7572 


4.8732 


133.5381 



